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TRANSIENT  TRANSPORT  IN  NANOSTRUCTURE 
DEVICES  VIA  THE  QUANTUM  LIOUVILLE 
EQUATION  IN  THE  COORDINATE  REPRESENTATION 


Abstract 

This  document  summarizes  work  performed  under  US  Navy,  Office  of  Naval  Research 
Contract:  N00014-95-C-0024  to  examine  transient  transport  in  quantum  structures  via  the  quan¬ 
tum  Liouville  equation  in  the  coordinate  representation.  For  comparison  transient  results  via  the 
Wigner  distribution  are  included.  The  issue  of  dissipation  in  nanostructures  was  also  initiated  in 
this  study.  A  summary  of  the  dissipation  results  is  included  in  a  reprint. 
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Introduction 


In  the  past  decade  there  have  been  demonstrations  that  band  structure  engineered  devices 
with  nanoscale  dimensions  have  the  potential  of  providing  a  dramatic  shift  in  semiconductor  tech¬ 
nology.  One  significant  illustration  of  this  is  the  development  of  multiple  value  logic  and  memory 
circuits  based  on  utilization  of  the  negative  differential  resistance  properties  of  either  isolated 
resonant  tunneling  structures  or  embedded  RTDs  (as  in  the  case  of  resonant  tunneling  bipolar  tran¬ 
sistors).  Another  is  the  demonstration  of  the  RTD  as  a  high  frequency  relaxation  oscillator,  sim¬ 
plifying  the  prospect  of  developing  high  frequency  clocks.  But,  the  operation  of  either  logic  or 
memory  devices  requires  knowledge  only  of  the  static  properties  of  RTDs,  as  does  the  operation  of 
the  RTD  as  a  relaxation  oscillator. 

But  as  device  size  shrinks,  the  active  region  of  the  device  approaches  the  size  of  the  elec¬ 
tron  wave-packet,  transit  times  approach  coherence  times,  and  the  device  operation  cannot  reliably 
be  treated  within  a  steady  state  picture.  Rather  transients  are  required  for  a  basic  understanding  of 
the  device  operating  principles.  This  study  under  US  Navy,  Office  of  Naval  Research  Contract: 
N00014-95-C-0024  was  concerned  with  transient  quantum  transport.  The  vehicle  for  examining 
transient  issues  was  the  density  matrix  in  the  coordinate  representation.  More  recent  studies  have 
focused  on  transient  transport  via  the  Wigner  distribution  function. 

Transient  transport  via  the  quantum  Liouville  equation  in  the  coordinate  representation  was 
begun  under  this  contract.  A  comparison  to  the  Wigner  studies  is  given  below.  Unless  otherwise 
noted  all  transient  calculations  are  coupled  to  Poisson's  equation. 

In  addition  a  study  of  dissipation  was  initiated.  This  is  summarized  in  a  reprint  included 
with  this  report. 

2.  Transient  Behavior  with  the  Time  Dependent  Liouville  Equation  in  the  Coordinate 
Representation  (Restricted  to  the  Single  Time  Representation) 

All  of  the  time  dependent  simulations  were  in  the  single  time  representation.  Two  time 
transients  were  formulated  but  not  implemented.  For  the  single  time  density  operator  the  quantum 
Liouville  equation  is: 

( i >  a> i&p- = O<o,  p„(o] + [h <»)***. , p„ <o] 

Here:  pop  (t)  is  the  density  operator  and  H{t)dissipation  is  that  portion  of  the  Hamiltonian  describing 
dissipation. 

Under  this  study  the  quantum  Liouville  equation  was  studied  within  the  coordinate  repre¬ 
sentation.  The  full  coordinate  representation  is  a  six  dimensional  space  plus  time,  where  the  den¬ 
sity  matrix  is  expressed  as:  p (x,  y,z‘,x\y',z';t).  We  dealt  only  with  the  restricted  situation  where 

the  carriers  are  free  in  two  directions,  y  and  z.  This  is  equivalent  in  the  Wigner  picture  to  dealing 
with  one-dimensional  transport  and  ignoring  any  non-parabolic  contributions  from  transverse  mo¬ 
mentum  states. 

Dissipation  was  introduced  through  quasi-Fermi  levels  as  shown  below  in  equation  (2): 


4 


dp(x,x',t)  _ 
dt 

p(x,  x\  t)  +  [(V(JC)  -  V  (x '))  -  (Ef  (x)  -  Ef  (x  ') )]  p(x,  X ',  t ) 

As  we  are  only  interested  in  differences  in  potential  energy,  we  are  only  interested  in  differences  in 
the  quasi-Fermi  levels.  Here: 

(3)  Ef(L)-Ef(0)  =  -\L  dxmv(x)r(x) 

In  the  above  the  velocity  v(jc),  is  obtained  as  the  ratio  of  the  position  dependent  particle  current 

density  and  the  carrier  density.  T(x)  is  a  position  dependent  scattering  rate.  The  particle  current 

density  is  given  by  the  diagonal  component  of  the  current-density-operator  in  the  coordinate  repre¬ 
sentation: 


ft2  (  <92  <92 

2 m  x2  d  x'2 


J(x,x\t ) 


2 mi  \  d  x  d 


-  pc*,* ',0(4) 


The  quasi-Fermi  energy  was  computed  subject  to  the  constraint  that  the  kinetic  energy  of 
entering  and  exiting  carriers  are  equal.  (The  quasi-Fermi  level  approach  to  dissipation  was  exam¬ 
ined  in  detail  in  a  variety  of  circumstances.  In  particular,  under  conditions  appropriate  to  classical 
transport,  all  of  the  standard  results  were  obtained.  In  particular.  Ohm’s  law  was  retrieved.  Under 
conditions  appropriate  to  quantum  corrections  this  dissipation  model  appears  to  provide  reasonable 
current- voltage  characteristics.) 

In  all  transient  calculations,  those  via  the  density  matrix  in  the  coordinate  representation 
and  those  via  the  Wigner  representation,  total  current  included  both  displacement  and  particle 
conduction  contributions;  the  former  proportional  to  the  time  derivative  of  the  local  electric  field. 
Thus  accurate  solutions  to  Poisson’s  equation  were  needed. 

The  boundary  conditions  on  equation  (2)  incorporated  the  density  matrix  equivalent  of  a 
displaced  Fermi-Dirac  distribution.  As  discussed  in  the  first  reprint,  the  upstream  boundary  condi¬ 
tion  is  p(x,x')cxp[i(x-x')mJ /(p0h)]  where p(x,x')  is  the  zero  current  quantum  distribution 

function  on  the  boundary,  p0  is  the  density  on  the  boundary,  and  J  is  the  current  density  at  the 

boundary.  This  condition  is  also  imposed  for  the  time  dependent  studies.  This  condition  is 
equivalent  to  imposing  displaced  Fermi-Dirac  conditions.  In  the  Wigner  function  calculations,  dis¬ 
placed  Fermi  conditions  were  not  imposed.  Rather  the  normal  spatial  derivative  was  set  to  zero. 
In  both  the  density  matrix  and  Wigner  function  calculations,  the  boundary  conditions  chosen  were 
designed  to  provide  flat  band  conditions. 


3.  Results  with  the  Density  Matrix  in  the  Coordinate  Representation 


For  the  problems  discussed  below  with  the  potential  energy  at  the  emitter  set  to  zero,  the 
system  is  brought  to  steady  state  for  a  given  collector  potential  energy.  A  step  change  to  a  new 
constant  value  is  introduced  for  V  (L) ,  and  the  resulting  time  dependent  behavior  is  computed  at 

fixed  increments  in  time.  All  calculations  assumed  Fermi  statistics,  parameters  appropriate  to 
GaAs,  except  for  the  barriers,  and  a  constant  effective  mass.  For  each  study  a  steady  state  was 
reached  at  a  value  of  100  meV.  A  step  change  in  voltage  to  150  meV  was  then  applied. 
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Several  structures  were  studied,  all  120nm  in  length.  This  device  length  is  too  small  for 
the  requirements  of  flat  band  at  the  boundaries  to  be  met.  The  choice  was  based  on  achieving 
small  time  steps  and  these  calculations  did  indicate  a  capability  to  perform  transient  studies.  The 
transient  studies  included: 

(1)  An  N+N"N+  triple  barrier  200  meV  structure; 

(2)  An  N+N'N+,  300  meV  double  barrier  structure,  and 

(3)  A  uniform  N-type  300  meV  double  barrier  structure. 

The  background  doping  and  structure  of  the  double  and  triple  barrier  results  are  shown  in 
figure  1 .  The  results  for  the  triple  barrier  structure  are  discussed  first. 


Figure  1:  Background  doping  distribution  and  the  double  and  triple  barrier  structures  for 
the  transient  calculations  with  the  density  matrix  in  the  coordinate  representation. 


The  triple  barrier  structure :  At  100  meV  the  steady  state  current  was 

1.477xl09 amps/m2.  Increasing  the  bias  to  150  meV  yielded  an  increased  steady  state  value  of 
current  equal  to  2.657x10 9  amps/ m2 .  Due  to  time  limitations  transitions  to  lower  current  with 
increased  bias  were  not  obtained  during  this  transient  study.  They  awaited  the  more  recent  Wigner 
calculation,  which  is  discussed  below.  The  time  dependent  space  charge  distribution  following  the 
step  from  100  meV  to  150  meV  was  followed  through  a  time  of  420  fs.  The  time  step  was  70fs, 
and  the  transient  current  at  the  end  of  each  time  increment  is  shown  in  figure  2.  While  we  did  not 
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carry  this  calculation  to  very  long  times,  the  initial  structure  in  the  calculation  shows  a  period  of 
near  200  fs.  This  number  is  significant,  as  we  see  very  similar  structure  in  transients  associated 
with  double  barrier  structures  obtained  with  the  Wigner  formulation. 

Examining  the  space  charge  distribution  following  the  step  change  in  voltage,  we  found  for 
the  cladding  regions  surrounding  the  barriers,  little  difference  in  density  at  the  different  bias  levels. 
Since  there  are  changes  in  voltage  in  the  surrounding  regions,  the  dominant  time  dependent  behav¬ 
ior  is  likely  due  to  displacement  current  contributions.  The  situation  is  different  in  the  interior. 

Figure  3  displays,  for  the  central  20  nm  region  the  initial  steady  state  solution,  the  final 
steady  state  solution,  and  for  two  early  time  intervals.  Several  points  are  noted. 

(1)  The  lowest  values  of  charge  density  occur  within  the  barriers,  and  there  are 
three  barriers.  (This  results  is  expected.) 

(2)  There  is  charge  accumulation  within  the  two  quantum  wells.  With  the  greatest 
amount  in  the  quantum  well  closest  to  the  emitter  contact. 

(3)  The  maximum  change  in  the  charge  distribution  within  the  two  quantum 
wells  appears  to  have  occurred  within  the  first  70  fs  time  step.  It  is  within  the 
quantum  wells  that  the  peak  carrier  density  occurs.  The  situation  within  the  bar¬ 
riers  is  different.  Here  there  is  considerable  transient  behavior  with  the  density 
not  quite  reaching  its  steady  state  value.  The  results  suggest  distributive  behav¬ 
ior  during  the  initial  transient.  Distributive  behavior  means  that  within  the  bar¬ 
riers  where  there  is  considerable  temporal  variation  in  the  charge  distribution, 
we  may  be  seeing  a  transient  inductive  contribution  to  the  potential  drop  in  the 
device.  In  the  quantum  wells  we  may  be  seeing  a  capacitive  contribution. 


0  70  140  210  280  350  420 


Time  (fs) 

Figure  2:  The  current  transient  for  the  triple  barrier  structure. 
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The  last  point  is  a  new  and  unexpected  result,  and  bears  careful  study  in  the  Wigner  formu¬ 
lation.  From  a  device  operational  point  of  view  this  means  that  most  of  the  longer  time  current 
transients  within  the  quantum  well  arise  from  displacement  current  contributions.  It  also  indicates 
that  there  is  considerable  physics  in  the  sub  70  fs  time  scale. 


Transient  Density  Following  Sudden  Bias  Change 


Distance  (nm) 


Figure  3:  The  transient  density  for  the  triple  barrier  structure  following  a  sudden 
change  in  bias. 


The  Double  Barrier  Structure :  The  situation  for  the  N+N"N+  double  barrier  structure  is 
similar.  Here  at  a  bias  of  100  meV,  the  steady  state  current  is  6.102  x  108  amps/m2,  at  150  meV 
the  current  is  1.047  x  109  amps/m2.  The  transient  density  distribution  shows  considerably  more 
structure  in  the  cladding  regions,  than  for  the  triple  barrier  structure  and  suggests  that  conduction 
current  contributions  as  well  as  displacement  current  contributions  are  significant.  The  situation  in 
the  interior  is  displayed  in  figure  4  and  is  very  similar  to  that  of  the  triple  barrier  structure.  Appli¬ 
cation  of  the  sudden  change  in  bias  results  in  an  increase  in  charge  everywhere  within  the  struc¬ 
ture.  But  here  the  transient  behavior  is  more  dramatic  than  that  arising  in  the  triple  barrier  struc- 
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ture.  First  we  see  similar  behavior  in  the  quantum  well,  where  the  carrier  density  appears  to  ap¬ 
proach  its  steady  state  value  in  a  time  less  than  70fs.  The  transient  behavior  associated  with  the 
barriers  is  much  different.  As  in  the  case  of  the  triple  barrier  structure  where  most  of  the  time  de¬ 
pendent  variation  in  charge  occurred  within  the  barriers,  here  the  time  variation  in  charge  density 
in  the  barriers  is  extreme.  The  consequence  of  the  barrier  charge  distribution  needs  to  be  exam¬ 
ined. 


Transient  Density  Following  Sudden  Bias  Change 


Figure  4:  The  transient  density  for  the  double  barrier  structure  following  a  sudden 
change  in  bias. 


The  Uniformly  Doped  Double  Barrier  Structure:  The  situation  for  the  double  barrier 
structure  embedded  in  a  uniformly  doped  1024/m3  region  also  shows  time  dependent  behavior,  but 
the  initial  and  final  state  distributions  of  charge  are  not  significantly  different,  and  the  steady  state 
distribution  of  charge  appears  to  be  reached  in  the  early  time  stages.  It  appears  that  dielectric  re¬ 
laxation  may  be  playing  a  significant  role  here. 

4.  Transients  Via  the  Wigner  Distribution-A  Comparison 
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The  detailed  charge  distribution  associated  with  barriers  is  far  richer  than  we  would  expect 
from  simple  time  independent  arguments.  More  recent  work  with  a  Wigner  transient  algorithm 
developed  at  SRA  indicates  some  unusual  features;  unusual  in  the  sense  that  they  are  unlikely  to 
predicted  from  a  simple  analysis  of  the  initial  and  final  states.  We  show  these  results  below. 

The  calculations  were  for  a  resonant  tunneling  diode  with  250  meV  barriers,  5  nm  wide, 
and  with  a  5  nm  wide  separation.  The  current  voltage  characteristics  demonstrated  a  current  drop 
back  with  no  dc  hysteresis.  The  device  was  subjected  to  a  transient  voltage  pulse  and  allowed  to 
settle  into  equilibrium.  The  voltage  pulse  started  from  approximately  230  meV,  a  value  below  the 
threshold  voltage  for  negative  differential  conductivity,  to  approximately  310  meV.  The  time  de¬ 
pendence  of  the  charge  distribution  is  shown  in  figure  5.  Figure  5  is  a  display  of  the  one¬ 
dimensional  distribution  of  charge  at  different  instants  of  time.  The  structure  is  200  nm  long. 

We  know  from  the  static  calculations  that  the  distribution  of  charge  prior  to  switching  to 
the  low  current  state  consists  of  considerable  charge  accumulation  in  the  quantum  well.  After  the 
transition  to  the  low  current  state  there  is  significant  loss  of  charge  within  the  quantum  well,  with  a 
large  increase  in  charge  on  the  emitter  side  of  the  barrier.  The  time  dependent  behavior  displayed 
in  figure  5  is  consistent  with  this  result.  But  there  is  more.  It  appears  that  the  transition  to  the 
high/low  charge  distribution  within  the  quantum  well  occurs  early  in  the  transient.  If  this  is  the 
case  what  is  the  origin  on  the  high  frequency  current  transients  that  commonly  accompany  the 
switching  transient?  (We  note  that  we  have  also  seen  these  transients  in  the  density  matrix  simula¬ 
tions.)  A  careful  look  at  the  simulations  indicates  that  the  charge  distribution  on  the  emitter  side 
of  the  barrier  is  undergoing  considerable  oscillation,  more  than  that  within  the  quantum  well. 
While  the  distribution  on  the  emitter  side  of  the  barriers  affects  the  charge  distribution  within  the 
quantum  well,  the  time  dependence  appears  to  be  dominated  by  the  charge  on  the  emitter  side  of 
the  barrier. 

The  time  dependence  of  the  charge  on  the  emitter  side  of  the  barrier  is  somewhat  obscured 
in  figure  5.  In  figure  6,  we  have  taken  out  the  transient  further  (in  time),  but  have  broken  up  the 
plot  into  two  sections:  an  early  and  late  time  transient.  There  is  considerably  more  quantum  well 
charge  during  the  initial  transient.  Note  the  large  charge  excursions  on  the  emitter  side  of  the  bar¬ 
rier.  (The  emitter  is  on  the  right  hand  side  of  both  figures.) 

The  conclusion  of  the  Wigner  transients  and  the  density  matrix  transients  is  that  the 
structure  of  the  transient  oscillation  will  also  depend  upon  the  design  of  the  regions  outside  of 
the  quantum  well. 
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Figure  5:  Wigner  function  calculation  of  the  transient  density  vs.  distance  vs.  time  for  a  resonant  tunneling 
diode  subject  to  a  sudden  change  in  bias. 


Density 


Figure  6:  As  in  figure  5,  with  two  exceptions.  Transient  is  carried  out  for  a  longer  period  of  time,  and  the  display  is  for 
the  initial  and  final  sections  of  the  transient. 


5.  Summary  and  Publications 

The  calculations  discussed  here  demonstrate  the  richness  of  the  transient  phenomena. 
There  is  tunneling  through  a  barrier  and  into  the  quantum  well,  coupled  to  displacement  current 
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Density  /( 10  e  +  24/m  e  +  3) 


effects,  and  contributions  arising  due  to  the  magnitude  of  the  density  within  the  structure.  We 
have  not  attempted  to  extract  a  tunneling  time  from  these  studies.  In  actual  device  studies  the  tun¬ 
neling  times  would  be  dressed  by  dielectric  contributions  and  the  variability  of  the  self-consistent 
field.  Furthermore,  large  signal  transients  would  dominate  indicating  that  tunneling  times  would 
depend  on  the  initial  and  final  states  of  the  system.  The  details  of  the  result  will  also  depend  upon 
the  scattering  parameters. 

The  transient  density  matrix  study  was  summarized  in  a  paper  that  is  included  with  this  re¬ 
port  (Reprint  1). 

As  indicated  in  the  introduction  there  was  also  a  new  formulation  of  dissipation.  This  is 
included  as  Reprint  2. 
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Reprint  1 


SELF-CONSISTENT  TIME  DEPENDENT  SOLUTIONS  TO  THE  QUANTUM 
LIOUVILLE  EQUATION  IN  THE  COORDINATE  REPRESENTATION:  APPLICATION 
TO  BARRIER  STRUCTURES:  Appearing  in  Hot  Carriers  in  Semiconductors,  K.  Hess,  J.  Le- 
burton,  U.  Ravaioli,  eds  Plenum  Press,  NY  Page  421  (1995). 
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SELF-CONSISTENT  TIME  DEPENDENT  SOLUTIONS  TO  THE  QUANTUM 
LIOUVILLE  EQUATION  IN  THE  COORDINATE  REPRESENTATION:  APPLICATION 
TO  BARRIER  STRUCTURES* 


H.  L.  Grubin1  ,  T.  R.  Govindan1 ,  and  D.  K.  Ferry2 

'Scientific  Research  Associates,  Inc.;  Glastonbury,  Connecticut,  USA 
2 Arizona  State  University;  Tempe,  AZ,  USA 


INTRODUCTION 


We  report  on  transient  accurate  self-consistent  solutions  of  the  quantum  Liouville  equation 
in  the  coordinate  representation: 


dp(x,x',t) 
d  t 


n2  f  d2 

2m  d  x2 
v 


d2 
d  x'2 


\ 

p(x,x\t)  +  [{V(x)-V(x'))-(EF(x)-EF(x')]]p(x,x\t) 

J 


where  EF(x)  and  EF(x')  represent  quasi-Fermi  levels  subject  to  the  constraint  that  the  kinetic  en¬ 
ergy  of  entering  and  exiting  carriers  are  equal  [Grubin  (1995)].  Dissipation  is  included  and  cur¬ 
rent,  J  Total  ,  contains  both  displacement  and  conduction  contributions.  Under  time  independent 
conditions  the  conduction  current,  j,  for  a  scattering  rate  T  satisfies  the  condition  that: 

Ef  (L)  -  EF  (0)  =  —j  [  dxmT(x)  /  p(x) .  The  boundary  conditions  incorporate  current  via  the  den- 

JO 

sity  matrix  equivalent  of  a  displaced  Fermi-Dirac  distribution  [  Grubin  et  al.  (1993)]. 

For  this  discussion,  V(x=0)=0,  and  for  a  given  downstream  potential  energy  V(L),  the  sys¬ 
tem  is  brought  to  steady  state.  A  step  change  to  a  new  constant  value  is  introduced  for  V(L),  and 
the  resulting  time  dependent  behavior  is  computed  at  fixed  increments  in  time. 

THE  RESULTS 

Three  symmetric  120  nm  length  structures  were  studied:  (i)  N+N'N+  [N+=1024/m3  (30  nm 
long),  N'=1022/m3  (30  nm  long)],  three  200  meV  barriers  [4  nm  barriers,  4  nm  wells];  (ii)  N+N"N+, 
two  300  meV  barriers  [5nm  barriers,  6  nm  well]  and  (iii)  as  (ii)  but  uniformly  doped  with  N 
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=1024/m3  .  We  assumed  Fermi  statistics,  parameters  appropriate  to  GaAs,  except  for  the  barriers, 
and  a  constant  effective  mass.  While  the  length  of  the  structure  is  too  small  for  any  real  device 
studies,  the  choice  achieved  small  time  steps.  A  step  change  in  voltage  to  150  meV  from  a  steady 
state  of  100  meV,  was  applied. 

The  triple  barrier  structure’.  At  100  meV  the  steady  state  current  was  1.477  x  109 
amps/m2 .  At  150  meV  steady  state  yielded  a  current  value  of  2.657  x  109  amps/m2.  The  time  de¬ 
pendent  space  charge  distribution  following  the  voltage  change  was  followed  through  420  fs.  The 
time  step  was  70fs,  and  the  transient  current  at  the  end  of  each  time  increment  is  shown  in  figure  1. 
70  fs  after  application  of  the  step  change  in  voltage  the  charge  is  generally  increased  everywhere 
within  the  structure.  (This  is  true  for  the  double  barrier  structure  as  well.)  For  the  regions  sur¬ 
rounding  the  barriers  we  find  little  time  variation  in  density  although  the  cladding  region  voltage 
undergoes  changes,  suggesting  that  the  dominant  time  dependent  behavior  in  the  cladding  region  is 
due  to  displacement  current  contributions.  For  the  barrier  region,  matters  are  different;  the  time 
dependent  changes  in  voltage  are  reduced.  At  70  fs  after  application  of  the  voltage  step  the  charge 
within  the  barriers  exceeds  that  of  the  lower  bias  steady  state;  see  figure  2,  which  displays  denstiy, 
within  the  central  20  nm  region.  At  140  fs  there  is  a  decrease  in  charge  within  the  three  barriers 
and  a  corresponding  increase  in  charge  between  the  barriers.  These  results  are  consistent  with  an 
interpretation  that  a  fraction  of  the  charge  present  at  70  fs  has  tunneled  through  the  barriers  at  140 
fs.  Excess  charge  in  the  wells  would  subsequently  tunnel  back  toward  the  respective  barriers  at  a 
later  time,  etc.,  until  a  steady  state  is  reached. 


Figure  1.  Current  transient  for  the  triple  barrier  structure. 
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Transient  Density  Following  Sudden  Bias  Change 


Figure  2.  Transient  space  charge  distribution  for  the  triple  barrier  struture. 


The  Double  Barrier  Structures :  The  situation  for  the  N+N'N+  double  barrier  structure  is 
similar.  Here  at  a  bias  of  100  meV,  the  steady  state  current  is  6.102  x  10  amps/m  ,  at  150  meV 
the  current  is  1.047  x  10  amps/m  .  For  the  regions  surrounding  the  barriers,  the  transient  den¬ 
sity  distribution  shows  considerably  more  structure  in  the  cladding  regions  than  for  the  triple  bar¬ 
rier  structure  and  suggests  that  conduction  current  contributions  as  well  as  displacement  current 
contributions  are  significant.  For  the  barrier  region ,  figure  3,  charge  appears  to  tunnel  to  the  cen¬ 
ter  well,  where  the  peak  density  exceeds  its  steady  state  value.  Tunneling  out  of  this  region  in¬ 
creases  the  charge  in  the  barrier.  This  tunneling  into  and  out  of  the  quantum  well  may  be  driving 
the  time  dependence  of  the  device. 

The  situation  for  the  double  barrier  structure  embedded  in  a  uniformly  doped  10  /m  re¬ 
gion  also  shows  time  dependent  behavior,  but  the  initial  and  final  state  disributions  of  charge  are 
not  significantly  different,  and  the  steady  state  distibution  of  charge  appears  to  be  reached  in  the 
early  time  stages.  Dielectric  relaxation  may  be  playing  a  significant  role. 
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FigureS.  Transient  space  charge  distribution  for  the  double  barrier  structure. 

SUMMARY  AND  REMARKS 

The  calculations  discussed  here  demonstrate  the  richness  of  the  transient  phenomena. 
There  is  tunneling  through  a  barrier  and  into  the  quantum  well,  coupled  to  displacement  current 
effects,  and  contributions  arising  due  to  the  magnitude  of  the  density  within  the  structure.  We 
have  not  attempted  to  extract  a  tunneling  time  from  these  studies;  in  actual  device  studies,  the  tun¬ 
neling  times  would  be  dressed  by  dielectric  contributions  and  the  variability  of  the  self  consistent 
field.  In  most  of  the  calculations  the  current  reaches  its  steady  state  value  at  approximately  500  fs, 
but  weak  oscillations  beyond  this  are  expected.  We  have  not  carried  this  further.  The  details  of 
the  result  will  be  depend  upon  the  scattering  parameters. 
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DISSIPATION  AND  QUANTUM  TRANSPORT  SIMULATIONS  IN  NANOSCALE 
DEVICES:  Appearing  in  Superlattices  and  Microstructures  20,  531  ( 1996). 
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A  key  simulation  issue  is  the  development  of  a  dissipation  formalism  for  the  time  dependent 
density  operator  equation,  that  is  amenable  to  numerical  methods  and  incorporates  the  role 
of  the  environment  on  state  renormalization  and  dissipation.  Using  standard  super-operator 
calculus,  projection  operators,  and  separating  the  system  of  interest  from  the  reservoir,  the 
relevant  operator  equation  is  derived.  In  addition,  the  role  of  the  reservoir  on  renormalizing 
the  energy  spectrum  is  discussed. 

©  1996  Academic  Press  Limited 
Key  words:  dissipation,  quantum  transport,  Liouville  operator. 


1.  Introduction 

Numerical  simulation  of  quantum  structures  has  been  a  major  element  of  device  physics  studies.  A  key  issue 
in  these  simulations  is  the  development  of  a  formulation  for  dissipation  within  the  framework  of  the  quantum 
Liouville  equation.  Such  a  formulation  is  considered  below  with  the  discussion  based  upon  Haken  [1]  for  a 
quantum  system  (e.g.  electrons),  which  may  be  far  from  equilibrium.  The  system  of  interest,  denoted  by  5, 
is  kept  far  from  equilibrium  by  coupling  it  to  other  systems  with  which  it  exchanges  energy.  The  discussion 
below  is  concerned  with  coupling  to  a  reservoir,  denoted  by  R ,  whose  detailed  properties  are  often  of  interest. 

The  quantum  Liouville  equation  for  the  total  system,  S  plus  R  is: 

it^-  =  [H,pT}  =  HpT,  (1) 

ot 

where  the  ‘carat’  over  the  operator,  designates  a  super-operator,  as  reviewed  in  Ref.  [2].  The  super-operator 
here,  also  called  the  Liouville  operator,  is  a  commutator-generating  super-operator,  whose  use  simplifies 
the  algebra  associated  with  separating  the  reservoir  and  system.  We  use  super-operator  algebra  to  obtain  a 
quantum  Liouville  equation  for  the  reduced  density  operator  for  the  system  S ,  defined  as  the  trace  over  the 
eigenstates  of  the  reservoir,  ps  =  Tr^  pr,  (with  a  similar  definition  for  the  reduced  density  operator  for  the 
reservoir,  pR  =  Tr(S)  Pr)- 

With  Hs,  Hr,  and  HSr  denoting,  respectively,  the  Hamiltonian  of  S,  R,  and  the  interaction,  the  Liouville 
equations  for  ps  =  Tr(fl)  pT  and  pR  =  Tr(S)  pT  are,  respectively: 

=  [Hs,  Ps ]  +  Tr w[HSr,  PtI  (2) 

Ot 

0749-6036/96/080531  +  04  $25.00/0  ©  1996  Academic  Press  Limited 


20 


532 


Superiattices  and  Microstructures ,  VoL  20,  No.  4,  1996 


M-Zr  =  [Hr,Pr]  +  Tr(S)  [//**,  fir).  (3) 

ot 

The  presence  in  eqs  (2)  and  (3)  of  the  trace  over  pj  means  that  we  do  not  have  a  prescription  for  obtaining  an 
equation  for  the  reduced  density  operator.  We  deal  with  this  below.  Our  approach,  which  is  based  upon  the 
projection  operators  of  Mori  [3],  allows  for  a  general  decomposition  of  the  system  and  reservoir  even  when 
they  are  significantly  entangled ,  as  is  often  the  case  for  far-from-equilibrium  situations. 


2.  Projection  operators  and  the  time  dependent  Liouville  equation 

Starting  from  equation  (1),  we  consider  the  projection  operator  P  and  its  complement  Q ,  such  that:  P+Q  = 
1.  The  operator  P  is  chosen  such  that:  Ppj  =  R  Tr(/f)  pr  =  Rps >  Tr(/?)  R  =  1,  PR  =  R,  where  R  is  chosen  to 
represent  the  uncoupled  distribution  of  the  reservoir,  e.g.  /?->/*  =  exp[— fiHR]/Tr{R){exp[—fHR]}.  Thus 
the  effect  of  the  operator  P  is  to  take  the  density  (or  other)  operator,  average  over  all  the  non-equilibrium 
coordinates  of  the  reservoir ,  and  generate  a  density  operator  that  is  a  product  of  the  equilibrium  reservoir 
operator  and  the  reduced  density  operator  of  the  system  [4].  The  dimension  of  the  subsequent  operator  is 
the  same  as  that  of  pr.  It  is  easy  to  demonstrate  that  P2  =  P>  and  that  P  is  a  proper  projection  operator. 
Then,  starting  from  eqn  (1),  and  after  some  super-operator  algebra,  the  time  dependent  Liouville  equation 
becomes: 


Here, 


in ^  =  HsPs  +  HrsPs  ~  ^Tr(«  HRS{t)  f  drU(t,  t)Hrs(t)Rps(t) 
ot  n  J  o 

Tr<«>  HRS(t)  f  dT(PU(t,  t)Hrs(t)RPs(t )) 
fi  Jo 

+Lxr(fi)  HRS(t)  f* dx(Q  J‘  dr,U(t,x’)PH(x,)V(t,,x)HRS(r)Rps(.T)y 


Hrs 

V(t,T) 

U(t,  t) 


=  Ttw(HrsR ), 

=  exp -ijf  dx'H(x')Q, 

=  exp  —  f  dx'H{x '). 
ft  Jr 


(4) 

(5) 

(6) 

(7) 


In  eqn  (4),  all  terms  are  expressed  as  functions  of  the  reduced  density  matrix  whose  time  dependence 
implicitly  includes  that  of  the  reservoir.  The  equilibrium  reservoir  coordinate  appears  to  first  order  in  the 
second  term,  and  at  least  to  second  order  in  the  remaining  terms  (on  the  right-hand  side  of  eqn  (4)).  The 
intra-collisional  field  effect,  for  example,  arises  from  the  presence  of  the  total  Hamiltonian  in  the  exponential 
of  the  operator  U{t,  r).  Below  we  introduce  approximations  to  place  the  above  results  in  perspective.  We  have 
also  assumed  that  at  t  —  0,  when  the  interaction  is  initiated,  the  density  operator  is  separable  into  orthogonal 
operators  for  the  system  and  the  reservoir.  In  this  case,  the  time  dependence  of  eqn  (4)  is  unaffected  by  the 
initial  condition. 

The  last  three  terms  on  the  right-hand  side  of  eqn  (4)  represent  the  effect  of  scattering  by  the  fluctuations 
in  the  reservoir.  These  terms  all  contribute  to  an  equivalent  self  energy.  On  the  other  hand,  the  second  term  on 
the  right-hand  side  is  a  non-dissipative  correction  to  the  system  that  arises  from  merely  opening  the  system. 
Traditionally,  this  term  leads  simply  to  a  modification  of  the  eigenvalue  spectrum  within  the  system.  However, 

the  nature  of  H  Rs  can  dramatically  affect  the  entire  structure  of  ps .  We  examine  a  test  case  in  the  next  section. 
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Fig.  1.  Plot  of  the  diagonal  elements  of  the  density  matrix,  i.e.  the  square  of  the  magnitude  of  the  wave  function,  for  a  square  quantum 
dot  subject  to  different  reservoir  couplings.  A,  The  coupling  is  through  tunneling  barriers.  B,  The  coupling  is  via  conducting  wave  guides 
supporting  two  modes.  Clearly  the  nature  of  the  coupling  creates  a  different  density  matrix  for  the  system. 


3.  Open  systems  near  equilibrium 

We  consider  a  square  quantum  dot  that  is  coupled  to  the  reservoirs  through  point  contacts  (these  are 
observable  in  figure  IB),  if  the  quantum  point  contacts  are  closed  down  to  create  tunneling  barriers  (the  weak 
coupling  limit),  then  the  properties  are  quite  uniform  across  the  dot.  In  figure  1  A,  the  diagonal  terms  of  the 
density  matrix  in  the  coordinate  representation,  As(jc,  y;  jcy)  =  |^(jc,  y)|2,  are  presented  for  the  tunneling 
barrier  case.  If  we  vary  an  applied  magnetic  field,  the  transmission  function  exhibits  a  series  of  sharp  tunneling 
peaks. 

In  contrast,  if  the  point  contacts  are  opened,  then  entering  particles  form  a  collimated  beam,  which  excites 
a  particular  set  of  eigenstates  required  to  reproduce  the  semi-classical  orbits  [5].  The  excitation  is  not  limited 
to  a  single  eigenstate,  but  excites  a  specific  set  of  modes  [6].  It  is  clear,  by  comparing  Fig.  1A  and  IB,  the 
actual  values  of  p$  are  dramatically  affected  by  the  specific  details  of  the  interaction  with  the  reservoir. 

Another  important  point  is  that  the  energy  level  shifts  and  coupling  that  occur  by  opening  the  system  to  the 
reservoir  are  distinctly  different  from  the  scattering  properties,  which  lead  to  level  broadening.  The  former 
arise  from  the  second  term  on  the  right  of  eqn  (4),  while  the  latter  arises  from  the  last  three  terms  on  the  right 
of  eqn  (4).  Scattering  from  the  reservoir  can  either  cause  broadening  of  the  levels,  or  actually  work  to  stabilize 
the  regular  orbits  [7],  It  might  be  expected  that  the  latter  case  would  actually  lead  to  narrowing  of  discrete 
energy  levels. 


4.  Scattering  in  non-equilibrium  systems 

More  generally,  what  is  done  with  eqn  (4)  depends  upon  the  problem  of  interest.  For  the  case  where  the 
Hrs  is  linear  in  the  momentum  [8],  the  dissipation  may  be  expressed  in  terms  of  a  quasi-Fermi  energy  model 
[9].  Ahn  [10]  treated  the  reservoir  as  stochastic,  performed  a  time  average  on  eqn  (4),  and  obtained  a  quantum 
kinetic  equation  for  interacting  electron  and  hole  pairs.  Krech  etaL[  11]  developed  a  master  equation  to  study 
macroscopic  quantum  tunneling  of  charge  in  ultra-small  single  electron  tunneling  double  junctions. 

We  shall  consider  some  of  the  more  general  features  of  this  equation.  The  third  term,  which  was  obtained 
without  the  use  of  perturbation  theory,  includes  the  intracollisional  field  effect.  If  the  off-diagonal  elements  of 
the  density  operator  in  the  momentum  representation  are  taken  as  small  compared  to  the  diagonal  elements, 
we  obtain  the  quantum  kinetic  equation  discussed  by  Ferry  [2],  as  well  as  that  of  Kreiger  and  Iafrate  [12].  This 
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third  term  incorporates  energy  conservation,  but  only  within  the  system  S .  We  note  that  the  results  of  Ref.  [12] 
were  obtained  through  perturbation  theory  with  the  interaction  Hamiltonian  as  the  perturbation  estimate.  The 
results  here  are  not  totally  dependent  upon  approximations,  and  suggest  a  broader  applicability  of  the  results 
of  Ref.  [12].  Higher-order  terms  corresponding  to  the  detailed  role  of  the  reservoir  are  readily  accessible. 
These  will  be  discussed  in  more  detail  elsewhere. 
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